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Abstract. We report the results of molecular dynamics simulations of an off-lattice 

protein model featuring a physical force-field and amino-acid sequence. We show 
that localized modes of nonlinear origin (discrete breathers) emerge naturally as 
continuations of a subset of high-frequency normal modes residing at specific sites 
dictated by the native fold. In the case of the small /3-barrel structure that we consider, 
localization occurs on the turns connecting the strands. At high energies, discrete 
breathers stabilize the structure by concentrating energy on few sites, while their 
collapse marks the onset of large-amplitude fluctuations of the protein. Furthermore, 
we show how breathers develop as energy-accumulating centres following perturbations 
even at distant locations, thus mediating efiicient and irreversible energy transfers. 
Remarkably, due to the presence of angular potentials, the breather induces a local 
static distortion of the native fold. Altogether, the combination of this two nonlinear 
efi'ects may provide a ready means for remotely controlling local conformational changes 
in proteins. 
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1. Introduction 

Biopolymers such as proteins and nucleic acids fold into complex three-dimensional 
structures, whose shape is strictly connected to their biological function [I]. The 
conformation of such molecules can change dynamically, in turn modulating the 
function: for example, activation or inactivation of enzymes relies on specific structural 
modifications occurring at specific locations [21 [Sj H]. Typically, such changes are 
driven by either mechanical forces or by converting chemical energy into conformational 
rearrangements and thus into mechanical work. 

Proteins in physiological conditions are immersed in a thermal bath and therefore 
exhibit random thermal fluctuations. However, the biological function of a given 
biopolymer is often closely related to a particular kind of motion, typically involving 
larges-scale vibrations [3 El [71 El [9] or the fluctuations of entire hinge-domain 

units [iniiiiittans]. 

If collective, low-frequency modes have been traditionally assumed to describe 
functional patterns, there is growing evidence that high-frequency vibrations contain 
information on protein stability [H]. Importantly, fast modes are strongly localized, 
due to the geometric heterogeneity of protein structures. Typically, such vibrations are 
localized at extremely stiff regions, such as hinges, which in turn assume a prominent role 
in regulating protein stability and function, as also suggested by experiments It is 
also known that active sites of enzymes have a marked tendency to be located within the 
stiffest segments of the structures [151 [IS] , which provides a further intriguing motivation 
to investigate the connection between localized vibrations and the peculiarities of protein 
folds. 

When nonlinear effects are considered, the connection between the dynamics of 
localized modes and the details of the scaffolds becomes more interesting. Recently, it 
has been shown that band-edge normal modes (NMs) can be continued to nonlinear 
localized vibrations of high energy termed Discrete Breathers (DBs) [15], also known 
as intrinsic localized modes [T7]. DBs are time-periodic, spatially localized solutions 
emerging generically in discrete networks of nonlinear oscillators [TBI HQ], that have 
been also observed experimentally in many systems [IHl [201 [211 [22]- 

While DBs have been widely studied in spatially homogeneous systems [I9], the 
role played by inhomogeneity on their properties remains largely unknown. However, 
the dynamical behaviours of DBs discovered so far in protein models are remarkable. 
In the context of the nonlinear network model (NNM), it has been shown that DBs 
emerge spontaneously, upon surface cooling, at few specific sites, invariably within the 
stiffest regions [23]. Moreover, they are also able to self-excite at a target site upon 
injecting some energy at a different location, thus mediating high-yield energy transfer 
events [2ll [251 [26] . In a more simplified model, it has been shown that DB excitation 
lowers the free-energy barrier associated with a given enzyme-catalyzed reaction [27] . 
thus confirming the role of protein dynamics in reaction-rate enhancement by enzymes 
highlighted by recent experiments [28] . 
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However, DBs in proteins have only been found and characterized so far in extremely 
simplified models, either taking into account the three-dimensional folds but with no 
heterogeneity in the force constants [I5] , or with slightly more elaborate potentials but 
imposing crudely simplified geometries [291 EOl IM] . Importantly, all studies performed 
so far notably lacked (i) realistic inter-particle potentials and (ii) an explicit account of 
the amino acid sequence. In this paper, we make a step forward by investigating DB 
excitation and their properties in a realistic off-lattice model of protein dynamics with 
coarse-grained realistic interaction potentials and a three-code amino-acid sequence. 

The paper is organized as follows. In the next section we introduce our system 
and describe the model, along with an account of our simulation and analysis protocols. 
In section [3} we describe the emergence of DBs as numerical continuation of NMs and 
provide a thorough characterization of their properties. Finally, we summarize our 
findings and outline possible future directions prompted by our results. 

2. Model and numerical methods 

The 3D off-lattice protein model studied in this paper is a modified version of 
the one initially introduced in Ref. [32] and successively generalized to include a 
harmonic interaction between next-neighbouring beads instead of rigid bonds [33] . 
This model has been studied to describe thermally-driven folding and unfolding [321 
[Ml l35l [331 [361 [371 [38l [39] and, more recently, to reproduce mechanical manipulation 
experiments [iQl HH [121 [131 [HI [15]. It consists of a chain of L point-like monomers 
mimicking the residues of a polypeptidic chain with an associated aminoacid sequence 
coded by a three-letter alphabet: hydrophobic (B), polar (P) and neutral (A^). The 
intramolecular potential consists of four terms: a stiff nearest-neighbour harmonic 
potential, Vh, intended to maintain the bond distance almost constant, a three-body 
interaction Va, which accounts for the energy associated with bond angles, a four-body 
interaction Vd corresponding to the dihedral angle potential, and a long-range Lennard- 
Jones (LJ) interaction, Vlj, acting on all pairs i, j such that |i — j| > 2, namely 

Vnir-i^i+i) = a{ri^i+i - r^Y (1) 
VA{e,) = A cos{e,) + B cos(2e,) - Vo (2) 

+ D,[l-S{e,,e,+i)cos{3^,))] (3) 

VlAt,,,) =^^A^-^] (4) 

Here, j is the distance between the z-th and the j-th monomer, while 6i and cpi are 
the bond and dihedral angles at the i-th monomer, respectively. The parameters a and 
ro = 1 fix the strength and the equilibrium distance between consecutive monomers 
along the backbone, respectively. The term VA^Oi) it such that it corresponds, up to the 
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second order, to a harmonic term kg{6i — 6q)'^/2, where 

_ cos Op _ ke 

'sin^^o' 4sin2^o' 
Vq = A cos Oo + B cos(26'o) , 

with ke = 20, Oq = bn/U rad or 75° gS]. 

The dihedral-angle potential is characterized by three minima for = (associated 
with the so-called trans state) and = ±27r/3 (corresponding to cis states). This term 
is mainly responsible for the formation of secondary structures. In particular, large 
values of the parameters Cj, Di favor the formation of the trans state and therefore of 
/3-sheets, while when cis states prevail a- helices are formed. The parameters (Cj, Di) 
have been chosen as in Ref. |36], i.e. if two or more beads among the four defining ip 
are neutral (N) then Cj = and = 0.2, in all the other cases Ci = D.^ = 1.2. The 
tapering function S{9i,9i+i) has been introduced in the expression of Vd in order to 
smooth off the singularities introduced by the derivatives of the dihedral angle cosines 
(for an exact definition see Refs. |l3l HZ]). 

The Lennard- Jones term Vlj describes the interactions with the solvent, that 
depend on the nature of the interacting residues as follows: if any of the two monomers is 
neutral, the potential is repulsive cn,x = and its energy scale is fixed by eN,x = 4; for 
interactions between hydrophobic residues cb^b = 1 and e^.B = 4; for any polar-polar 
or polar-hydrophobic interaction cp^p = cp^p = —1 and ep^p = ep^p = 8/3. 

According to the above definitions, the Hamiltonian of the system reads 

i=i ^ i=\ 

+ E ^a(^.) + E ^d(<^., + E E yLJ{r^,^ (5) 

i=1 i=2 1=1 j=i+3 

where all monomers are assumed to have the same unitary mass. 

In the present paper we consider the following sequence of L = 46 monomers: 
BgN^[PB)/iN:iBgN^[PB)^P . This sequence has been widely analyzed in the past for 
thermal folding [321 EH ESI ESI E31 EZl ESI EH] as well as for mechanically induced 
unfolding and refolding [iQl HH HJl |l3]. Here, we adopt the same potential and 
parameters as in Refs. [Ml 1121 113], except for a stiffer harmonic constant a (50 < 
a < 1000). The simulations reported henceforth refer to a = 1000, except when 
otherwise indicated. We have verified that this choice does not affect the thermodynamic 
properties of the model {i.e. the folding temperature and the hydrophobic collapse 
temperature) found in Ref. [13] for a = 50. This result could be expected, since the 
harmonic bead-bead interaction introduced in [33] has been already shown not to alter 
the folding properties of the original Honeycutt-Thirumalai model [32] . 

With the above choice of parameters, the heteropolymer exhibits a four-stranded 
/3-barrel native configuration (NC), described by the coordinates qNc{i), i = I ■ ■ ■ L, and 
shown in Fig. [1} The NC corresponds to the absolute minimum of the potential energy, 
whose value for a = 1000 is VJvc ~ —48.94. The native structure is stabilized by the 
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Figure 1. (Color online) Native configuration of the protein model. Different colors 
represent the three different types of beads: white (neutral), light blue (polar) and 
sand (hydrophobic). Parameters are: L — 46, a — 1000. 

attractive hydrophobic interactions among the B residues, in particular the first and 
third Bg strands, forming the core of the NC, are parallel to each other and anti-parallel 
to the second and fourth strand {PB)^ and {PB)^P. The latter are exposed to the 
solvent due to the presence of polar residues. Overall, the four strands are separated by 
stretches of three consecutive neutral beads, forming as many turns (see Fig.|l|). These 
involve the following sites: 10-12 (first turn), 21-23 (second turn) and 33-35 (third turn). 

In the following, we report equilibrium microcanonical results obtained by 
integrating Hamilton's equations by means of a fourth order symplectic integrator |18] 
with a time-step of 10~^ time units (ensuring a relative energy conservation of ~ 10~^). 
In all cases, the initial coordinates of the beads were taken to correspond to the NC 
(zero-displacement) . 

2.1. Normal modes analysis 

The properties of nonlinear modes in spatially and force-heterogeneous systems strongly 
depend on the features of the linear spectrum: in particular, non-resonant nonlinear 
modes may emerge within inter-mode gaps and highly localized vibrations can be 
associated with band edges [15]. Hence, a detailed understanding of the hnear spectrum 
is an essential pre-requisite for our analysis. More precisely, we wish to investigate how 
the normal mode frequencies and eigenvectors depend on the harmonic bond stiffness 
a. 

The NM spectrum associated with the NC of our model protein is composed by 
3L — 6 non-zero frequencies {fk = A/Afc/27r}, A^, being the eigenvalues of the Hessian 
matrix |19]. The NM spectra are reported in Fig. [2] (a) for values of a ranging from 50 to 
1000. The most important observation for our purposes is that, for large enough values 
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of a, the spectrum splits in two bands separated by a gap. The highest portion (the 
optical band) comprises L — 1 modes associated with the stiffest force constants, that is 
bond (backbone) distortions along the chain, while the lowest part (the acoustic band) is 
composed by the remaining 2L — 5 modes that represent collective motions of the beads 
which are hardly affected by changes in a. This interpretation can be substantiated by 
the following analysis. We have first slightly perturbed the NC along the direction of 
each eigenvector with an arbitrary amplitude e, i.e. qNcij) ^fe(^) = QNcii) + fc(0- 
Then we have evaluated the relative variations of the potential energy components Vj, 
j = H (harmonic), A (three-body angular), D (diehdral), LJ (Lennard- Jones), with 
respect to their values in the NC according to the following definition: 

. \Vjiqk) - VjiqNc)\ 

As it is clear from Fig. [3} the perturbation along the first 45 eigenvectors is only restricted 
to the degrees of freedom associated with bond deformation, i.e. to the harmonic 
contribution to the potential energy. 
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Figure 2. (Color Online) (a) Frequencies of the NMs, fk- (b),(c),(d),(e) Inverse 
of the participation ratio for the corresponding eigenvectors for increasing value of 
the stiffness of the harmonic bonds. Triangles, squares, stars, and circles refer to 
a = 50, 100, 250, and 1000, respectively. 



Furthermore, the formation of the gap is accompanied by a localization of the 
eigenvectors around the gap- facing band edges. In order to characterize the degree of 
localization of the k-th NM e ^(i) (i = 1, 2, . . . , L), we measured its inverse participation 
ratio [50j, 

^k = j:\e,{^t (7) 

■t=i 

where the eigenvectors are normalized to unity. For an eigenvector localized on a single 
site ^ ~ 1, while for a completely delocalized state ^ ~ 1/L. Therefore, the more the 
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Figure 3. (Color online) Relative variations of the energetic contributions due to a 
perturbation of the NC along the direction of the fc-th NM, as given by eq. for 
a = 1000. Circles, triangles, diamond and stars correspond to Afc^Vtr, A^^Va, A^^Vd, 
Afe^y^ ,, respectively. 

eigenvector is localized the larger is ^. The inverse participation ratio is plotted in Fig. |2] 
(b-e) for different values of a. In particular, for a = 50 all the eigenvectors are extended 
(see Fig. [2] (b)). However, by increasing a, the degree of localization of NMs in the lower 
optical band edge {k = 43, 44, 45) and of the first three NMs of the acoustic band edge 
(k = 46, 47, 49) increases (Fig. |2] (c-e)). A further increase in stiffness does not only lead 
to an enhancement of localization in the proximity of the gap, but also determines the 
localization of the first three optical eigenvectors {k = 1,2,3), see Fig. [2] (e). It is also 
interesting to observe that for increasing a the first and last three optical frequencies 
detach more and more from the band core (see again Fig. [2] (a)). In the following, we 
will fix a = 1000. 

We notice that the most localized modes of the optical band are localized at the 
three turns of the NC. More precisely, modes = 1,45 at the first turn (Fig. |4]), NMs 
= 2,44 at the third turn and modes /c = 3,43 at the second one. The modes at both 
edges of the optical band are characterized by single or groups of neighbouring oscillators 
in opposition of phase, as shown in Fig. |4j Qualitatively, one may interpret them as 
impurity modes with the turns acting as structural defects [51]. Interestingly, for a 
toy model reproducing a single bent 2D chain with fixed curvature, Archilla et al |30] 
also reported a pair of localized modes lying at the band edge, with frequencies slightly 
detached from the band itself. Concerning the acoustic band-edge modes k = 46, 47, 48, 
only the y-component appears localized on the first and third /3-strand Bg, with the 
oscillators on the first and third strands in opposition of phase. Since these strands 
represent the core of the protein, we expect that excitation of these modes should cause 
large rearrangements in the structure, leading to protein unfolding. 




Figure 4. (Color Online) Cartesian components of the eigenvectors e i (top) and e 45 
(bottom) corresponding to the frequencies at the upper and lower edges of the optical 
band, respectively, for a = 1000. 



3. Emergence of discrete breathers 

In this Section we show how DBs can be created by exciting the NC along the direction of 
certain NMs. In these simulations, we initialize the beads positions in the NC and assign 
the initial velocities proportional to the pattern of the selected NM. The amplitudes of 
the kinetic energy perturbation Kq will be henceforth measured in temperature units 
T = 2Kq/3L, where a unitary Boltzmann constant is assumed. 

As we shall show, the excitation of a normal mode provides an effective means of 
feeding energy to a discrete breather. However, due to specific selection rules matching 
spatial overlap between normal modes [M], a portion of the initial energy necessarily 
fiows into a number of other modes. Such background radiation competes with the 
nonlinear localization process. This in turn, as we shall see in the following, accelerates 
the collapse of nonlinear modes. In order to analyze DB properties, it proves useful to 
get rid of background radiation through surface damping To this end, in the first 
part of the work we shall cool down the velocities of the beads at the chain terminals. 
More precisely, we integrate Hamilton's equations of motion while rescaling, at regular 
time intervals At, the momenta of the beads located at the chain terminals i = 1, . . . , 4 
and z = 43, . . . , 46), that is 

p{i) -> p '{i) = (1 - l)p{i) with 7 < 1 . (8) 

Unless otherwise specified, the cooling parameters are fixed to 7 = 0.001 and At = 0.02. 
In section 3.1.1| a comparison between the dynamics with and without cooling is 



presented. 

In order to visualize the energy localization along the chain highlighting the 
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Figure 5. (Color Online) (a) Total energy difference per site AEtotii) with i = 
1, . . . ,L, as a function of time, for DBl with an initial temperature T — 0.63. (b) 
Kinetic site energy (i?fci„(i)) averaged over a time span of the order of 50 — 100 
breather's periods in order to smooth out fluctuations. From top to bottom DBl 
with initial temperature T = 0.63, DB2 at T = 0.42 and DBS at T = 0.63. Cooling 
parameters are 7 = 0.001, At = 0.02 (a) and At — 0.005 (b). Grey regions mark the 
protein turns (see text). 



emergence of a DB, we recorded the site kinetic energy Ekin{i) and the total excitation 
energy per site AEtot{i) = Ekinii) + V{i) — Vi^cif), where V{i) is the potential energy 
contribution which can be associated with the i-ih bead (see Appendix A for the exact 
definition) and VJvc(^) is the corresponding value in the NC. The presence of a DB is 
detected as a protracted concentration of energy on a limited number of sites. 

Our first important observation is that a stable DB can be created by perturbing 
the NC along the direction of each of the first 45 modes, or along a linear combination of 
them. An example is reported in Fig. |5] (a) for the excitation of the lower optical-band 
edge mode. Remarkably, only three distinct breathers were observed to emerge despite 
the different directions of the employed perturbations. They are located close to the 
three turns and we thus denote them as DBl, DB2, and DBS, after the index of the 
corresponding turn (Fig. [s] (b)). 

A peculiarity of such localized modes is that they feature patterns that are 
not simply localized over a few adjacent sites. Rather, they display non-negligible 
components localized elsewhere in the chain. This can be appreciated by looking to the 
average site kinetic energy reported in Fig. |5] (b) for the three breathers. The largest 
contribution to the DB total energy comes indeed from the few sites located at the 
corresponding turn. However, energy components two-orders of magnitude smaller can 
be found over the other turns too. The origin of this peculiar localization patterns 
can be found in the three-dimensional structure of the protein, and it is induced by the 
requirement of momentum conservation. In other words, the special built-in momentum- 
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Figure 6. (Color online) DB obtained as a continuation of NM k =45. (a) Time 
evolution of the projection 7145 for T = 0.21. The inset shows a close-up of the short- 
time dynamics, (b) 1 — 7f45 as a function of the total energy E^b- The dashed line is 
a power law of the type 1 — 7f45 oc E^b with 77 « 0.88. 



conservation properties of NM patterns are reproduced by DBs. As a matter of fact, 
it turns out that a mode that is locahzed at a few sites conserves its momentum by 
fractioning an equal and opposite amount of momentum among few other locations, 
instead of spreading it over the whole structure. The interesting question arises whether 
this is a generic feature of NMs in proteins or it is rather the expression of the peculiar 
native fold considered here. 

The three DBs can be considered as nonlinear continuations of three corresponding 
NMs, that are localized at the same turns and whose vibration patterns can be 
considered as precursors of the DB displacement fields. This is confirmed by looking 
at the normalized projections vr^ of the breathers' velocity fields on the corresponding 
modes, defined as: 




where the angular brackets denote the average over a window of w time units [w = 4.0 
in the following). Fig. [6] (a) shows that, for DBs originating from NM 45, after an initial 
transient 7145 shows a stable trend with very small residual oscillations around the mean 
value 7f45 evaluated in the last stage of the simulation (3000 < time < 4000). The 
projections on the other modes (not shown in the figure) are smaller than 0.015. In 
Fig. [6] (b) we plot 1 — 7145 as a function of the DB total energy, E^^g = K + V — VJvo 
measured in the last stage. The figure shows that 1 — 7145 vanishes when the energy is 
decreased following a nearly linear trend. 

These results confirm and generalize what reported in the context of the Nonlinear 
Network Model (NNM), where gap-less breathers were shown to arise as continuations 
of edge NMs [15]. Also in the framework of the NNM, few special regions were shown 
to act as energy-accumulating centres upon generic excitation of the system, exactly as 
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we observed here [23] . 

Perturbations along all the optical modes of amplitude T = 0.63 resulted in the 
excitation of DBl or DBS, while DB2 was observed only in a few cases. The reason 
why DB2 seems somehow more difficult to excite is likely to depend upon the peculiar 
structural neighborhood of the second turn. This region lies deeper inside the core of 
the protein, where beads are more constrained, thus hindering the DB oscillations (see 
also the cartoon showing DB patterns in Fig. Isl) . 




^DB 



Figure 7. (Color online) DBs originating from k =45 (DBl), k =43 (DB2) and k =44 
(DBS). Frequencies of DBs as a function of their total energy Edb- The horizontal 
dashed lines flag the frequencies of the corresponding NMs. The inset displays E^b as 
a function of the initial energy Eq = Kq — Vjvc for the three DBs: circles, triangles, 
squares refer to DBl, DB2, DBS, respectively. The dashed line marks the full-efficiency 
regime Eos = Eq. 



The DB frequencies fnBi^ fDB2 and /dbs are computed from the power spectrum 
of the displacements of the most energetic bead within each DB. Typical DB spectra 
feature a main peak and minor (several orders of magnitude smaller) peaks signalling 
the residual excitation of few normal modes. The measured frequencies are shown in 
Fig. [7] as functions of the DB total energy Edb, all three DBs exhibit an almost linear 
decrease characterized by the same dispersion relation fj^s^Ej^g) = fk{l — Ejjg/e), where 
fk is the frequency of the NM from which the DB originates and e ~ 250. We notice 
that e could be in principle calculated through Lindstedt-Poincare perturbation theory 
as done, for example, in Ref. [T5] . 

The decrease of the DB frequency with the energy Edb suggests that the relevant 
nonlinearity is of the soft type [IH] . This explains why the presence of a gap is necessary 
for the formation of localized nonlinear modes. Indeed, we were unable to excite DBs 
in cases where the spectrum is gap-less. 

The frequencies foBi and foBs originate from the bottom of the optical band 
(/44 ~ /45) and enter the gap for arbitrary small energies. On the contrary, /db2 
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lies between and /44 for small energies, staying manifestly non-resonant with the 
neighbouring optical modes. This behaviour is reminiscent of intra-band DBs predicted 
by Kopidakis and Aubry in nonlinear disordered chains [521 153] and found analytically 
in nonlinear network models of proteins [15]. For larger energies, the frequency decreases 
further and also these DBs eventually enter the gap. 

It is instructive to consider whether a given DB can be created with arbitrary 
energy by exciting a given normal mode. The inset in Fig. [7] shows the relation between 
the initial total energy Eq = Kq — VJvc and the energy that is found stored in the DB. 
The ratio of the two values can be interpreted as an estimate of transfer efficiency, i. 
e. the fraction of the initial excitation energy channeled into the DB and consequently 
pinned down at a very specific location. Remarkably, we see that all the three DBs are 
characterized by a nearly unitary efficiency, up to a point where the curve saturates and 
then starts bending down. The transfer of energy from the NM to the DB is nearly 
optimal up to initial energies of the order Eq = Kq — V^^ ^ 30. This indicates the 
existence of a sort of maximal "DB capacity" - if the NM is feed with larger energies, 
the excess energy cannot be further injected in a localized mode and it is necessarily 
spread across the whole structure. 

Finally, no DB could be excited by perturbing the structure along acoustic NMs. 
This is most likely a consequence of the soft nonlinearity. Moreover, we were also 
unable to excite DBs through arbitrary localized perturbations of the NC. Indeed, initial 
conditions of the latter type have non-zero components over many NMs, including the 
acoustic ones. 

3.1. DBs originating from the lower optical-band edge 

In this section we shall present a more detailed analysis of the properties of DBl orig- 
inating from the bottom optical mode k = 45. This breather is localized on the first 
turn, from site 10 to site 14, as it is shown in Fig. |5j The energy redistribution dynamics 
within the DB can be appreciated from Fig. [9] (a). In particular, it is clear that the 
two site pairs (12,13) and (11,14) exhibit approximately phase-locked oscillations, while 
they are anti-phase locked among them. This is another illustration of the fact that the 
DB exchanges little or no energy with the surrounding. 

It is instructive to analyse the different contributions to the total energies when 
the DB is present, in order to quantify to what extent the different degrees of freedom 
participate to the breather dynamics. As it can be seen from Fig. |9] (b), at one of the 
most energetic sites, i = 12, the DB vibration involves essentially the backbone bonds 
(^ 84 — 86 %) with a small contribution coming from the angular terms (~ 13 — 15%), 
while the dihedral and Lennard- Jones contributions are negligible. The observed ratios 
between the different energetic contributions do not change by varying the initial energy. 
We observe the same behaviour for all the sites in the interval (11 — 13). 
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Figure 8. (Color online) Displacement field of DBl. Side (a) and front (b) views 
(only displacements larger than 1.0 have been shown as arrows). The initial condition 
corresponds to T = 0.63. 
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Figure 9. (Color online) (a) Total site energies for DBl and (b) the different energy 
contributions at site z = 12 versus time {AEj = Vj — V«c(i) where j = H,A, the 
dihedral and Lennard- Jones contributions being negligible). Parameters are: T ~ 0.63, 
7 = 0.001, At = 0.005. 



A closer inspection of Fig. [9] (b) shows that the harmonic bond component oscillates 
at frequencies 2f£)B but also has a sizeable component at /db- This is due to the soft 
nonlinearity which is known to induce a DC component in the displacement patterns of 
localized modes |5H [19] . To understand this in a simple way, let us consider two bond 
lengths along the backbone at site n , ri = r„_i „ and r2 = r„ „+i. We may write their 
evolution as an oscillating part plus a static distortion of the bond lengths, 

ri(t) = Ai sm{ujt) + (ri) 
r2(t) = A2sm{ujt) + (rs) 
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where oo = 27rfDB and the angular brackets represent a time average. Using the 
definition reported in Appendix A, it can be shown that the contribution of the two 
bonds to the local harmonic energy is 

Vu{n)<x^-^i^^sm\ujt) 

+ [A((ri)-ro) + A2((r2)-ro)]sinM (10) 
{{n) - rpf + {{r2) - rpf 
2 

Thus, the presence of a component at frequency u implies that (ri) and/or (r2) must be 
different from their equilibrium value tq. Accordingly, we conclude that the emergence of 
a DB also causes a stable structural distortion of the protein. Given that Ai ^ A2 ^ 0.1 
(see again Fig.jo]), also tiny differences {ri)—rQ ^ 10~^, as we recorded in our simulations, 
are able to affect the local bond energies. 




Figure 10. (Color online) Difference between the average values of the bond-bond 
angles 6i when the DBl is present and the corresponding equilibrium values (in the 
NC), 9i,is!C: when the total energy of DBl is E^b = 34.6. The inset shows a close- 
up of the kink region for increasing energies: (red) squares, (blue) triangles, (green) 
diamonds, (black) circles refer to Eob = 26.3,31.68,33.61,34.6, respectively. 



Much more prominent is the distortion effect caused by the soft-nonlinearity in 
the angular degrees of freedom. In particular. Fig. 10 reveals that the breather is also 
characterized by a kink-shaped angular distortion. The kink's amplitude increases with 
the total energy stored in the breather (see the inset in Fig. 10). 



3.1.1. DB collapse and the effect of cooling In the case where no cooling is applied 
(7 = 0) the DBs have a finite lifetime. To illustrate this, in Fig. 11 (a) we report the 
time evolution of the projection 7145 of the DBl, obtained as a continuation of the lower 
optical-band edge mode {k = 45). The sudden drop of 7145 signals the collapse of the 
excitation, accompanied by a rapid redistribution of energy over all the normal modes 
(equipartition). The typical collapse time increases upon decreasing the initial energy 
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(compare the three curves in Fig. 11 (a)). This is similar to what demonstrated for 



one-dimensional chains where the time to equipartition scales as an inverse power of the 
energy density [55l [56] . 

The effect of the collapse of a DB on the protein structure can be appreciated by 
considering the Kabsch distance 6k, which is a commonly used measure of the structural 
distance between two protein configurations [57]. In particular, Fig. 11 (b) shows how 



the distance 6k of the protein structure from the NC evolves in time in the presence of 
a discrete breather of type DBl. As long as the DB is present, 6k fluctuates around a 
relatively small value {6k — 0.1), meaning that the protein structure does not deviate 
appreciably from the NC. In fact, the static distortion effect illustrated above only 
concerns a small region of the fold. On the other hand, when the DB collapses, 6k 
starts to fluctuate wildly around a substantially larger value {6k ~ 0.35), thus signaling 
the occurrence of important conformational rearrangements. 

The DB lifetime gets substantially increased by removing the background vibrations 
through cooling. Fig. 



11 



(c) clearly illustrates how dissipation renders the DB stable by 
quickly eliminating background radiation. The stabilizing role of the boundary cooling 
method is well documented [5H1 EHl [23] and our results confirm that it can be used 
conveniently also for such a complex structure. 
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Figure 11. (Color online) (a) Projection of the DBl velocity field on the fc = 45 
NM, 7r45, as a function of time for increasing values of the initial energy. From the 
top to the bottom curve, T = 0.42,0.63,0.84 (non-dissipative dynamics), (b) Kabsch 
distance Sk for T = 0.63 (non-dissipative dynamics). The vertical arrows mark the 
time at which the DB collapses, (c) Comparison of 7145 for dissipative (upper red line), 
and non-dissipative (lower black line) dynamics for T = 0.63. 
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3.2. Excitation of optical band-edge NMs: long-range energy transfer 



The excitation of NMs at the upper edge of the spectrum gives rise to quite pecuhar 



phenomena. As a first example, Figs. 12 depict the emergence of a DB obtained by 
exciting the first two highest-frequency NMs, i.e. Ci, which is locahzed at sites 10 — 12 
and 62, which is locahzed at sites 34 — 36. In the first stage, the energy quickly spreads 
and remains confined in the vicinity of the perturbed location. After some time, a DB 
self-localizes, harvesting energy from the background and pinning it in the same region. 
The DBs obtained are again DBl and DB3. This is reasonable in view of the similar 
spatial structure of Ci, 6*2 and DBl, DB3, respectively (see again Figjl]). Moreover, we 
found that the frequencies of DBs originating from ci lie on the same frequency-energy 
curve as the breather originating from normal mode 45, shown in Fig. [7j 
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Figure 12. (Color online) Time evolution of total site energies starting from excitation 
of the upper band-edge optical modes: (a) Emergence of DBl from ei, T = 0.63 and 
(b) of DB3 from es, T = 0.84. 



Remarkably, we have also observed that the excitation of an edge mode may result 
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in a long-range energy transfer event. An illustration of this phenomenon is given in 



Fig. 13 a), for the excitation of the third highest-frequency NM. The NM 63 is localized 
at sites 22 — 27. Immediately after the perturbation, the DB appears localized in 
the same region, but seems to collapse rapidly afterwards, spreading its energy evenly 
across the structure. However, after a considerable time span, another DB emerges at 
a different location, namely at sites 11 — 13. This is the region where the lower optical- 
band edge mode {k = 45) is localized. A substantial part of the initial energy has been 
transferred and pinned down irreversibly at the other end of the structure, covering the 
distance from a turn to the following one (see again the cartoon depicting the protein 
structure). 

Such energy transfer phenomenon can be rationalized by analyzing the projections 
of the DB velocity field on the NMs during the time evolution of the perturbation. From 



Fig. 13 (b) one can clearly see that the initially excited mode is quickly emptied of its 
initial energy, which gradually flows into other two modes. As it shows, asymptotically 
the two target modes almost describe the entirety (7144 + 7145 ^ 1) of the DB, which 
is in fact localized at the turn opposite to the excited one, where the two target NMs 
are also localized. The reverse process is forbidden since the DB frequency is no longer 
resonating with the original mode. We remark that this phenomenon, although evocative 
of resonant energy transfer among few selected normal modes in proteins, is a one- 
way transfer, as the DB will retain the transferred energy for times comparable to its 
lifetime [601 ED. 
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Figure 13. (Color online) Long-range energy transfer starting from 63 at T = 0.63. 
(a) Time-evolution of total site energies; (b) projections of the DB velocity field on 
three NMs as functions of the time. 



3.3. Excitation of deep-band NMs 

NMs in the core of the optical band are generally characterized by a low level of 
localization, as one can see in Fig. |2] (b-e). As a result, when one of such modes is 
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excited, we observe a complex localization pattern, that alternates site-hopping and 
energy-pinning stages to phases where the energy is more evenly distributed. After a 
certain time, a stable DB emerges, focusing the energy at one of the turns. A typical 
realization of this scenario is illustrated in Fig. [I4|(a) for the excitation of the fifth 
highest-frequency mode. 

A different phenomenology is observed when we excite a mode lying deeper within 
the band. It may occur that the excitation energy remains trapped in two DBs localized 
on two different turns, realizing a state which is reminiscent of multibreather states 
observed in non-linear disordered systems [521 [53] . A realization of this phenomenon for 



NM 33 is reported in Fig. 14 (b). In particular, we notice that starting from NM 33 and 
changing the initial energy we can get different solutions. For example, for T = 0.21 we 
obtain a solution localized on the first and second turn, while the multibreather localizes 
on the first and third ones for T = 0.63. 

A very peculiar localized solution develops through the excitation of normal mode 
k = 36, whose pattern is localized in the protein's tail. As it is shown in Fig. 



14 



the energy remains almost entirely confined to the edge sites 44 — 46, with smaller. 



but non-negligible, amounts of energy also involving the turns Fig. 15 (a). We term 



such excitations boundary breathers (BE). As in the case of DBl, DB2 and DB3, the 
frequency of oscillation of the BB also lies in the gap {/bb = 7.739 for E^^g = 34.75). At 
variance with DBl, DB2, DB3, the energy of BBs has not only harmonic and angular 
contributions but also a sizeable Lennard- Jones one, resulting from a non-negligeable 
interaction between the first and last beads, as well as the second turn region of the 
chain which are relatively close in the NC. We note that nonlinear surface states of 
different kinds are known and studied in many contexts [121 



Eventually, the BBs hop to one of the three breathers described above. Fig. 16 



illustrates an example of this process, where the BB transfers its energy to DB3, as 
confirmed by the time evolution of the normal-mode projections shown in Fig. 15 (b). 
Asymptotically, we see that the main projection is vrae, while at t ~ 2.75 x 10^ an abrupt 
transition occurs and the leading projection becomes 7144, signalling the energy transfer 
to DB3 (corresponding to NM k = 44) localized at the third turn. 

While at high energies, i.e. T = 0.84, 1.05, the lifetime of BBs before a transition 
to DB3 takes place is short, at lower energies we observe the formation of a transient 
multibreather-like state with a progressive transfer of energy from BB to DB3. One of 
such examples is depicted in Fig. 16 Decreasing the excitation energy further, BBs are 
no longer observable. 

We note that the end of the protein can be considered as a defect /discontinuity 
in the chain, since it is a free end that interrupts the chain itself. Therefore, BBs 
can be considered as akin to localized modes emerging in other contexts close to point 
defects. Furthermore, we were unable to excite stable BBs locahzed on the first beads 
of the sequence, despite NMs with significative components on that terminal do exist. 
A possible explanation of this apparently contradictory result is the fact that the first 
strand is part of the protein core. Therefore, it is a much more rigid structure with 
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Figure 14. (Color online) Evolution of total site energies starting from excitation 
of the deep-band optical modes: (a) Emergence of DBl from e^, T — 0.63 and 
(b) multibreather solution localized on the first and third turn starting from 633 at 
T = 0.63; (c) Boundary Breather (BB) following the excitation of mode 636, T = 0.84 
(non-dissipative dynamics). 
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respect to the last strand, which instead is exposed to the solvent and can oscillate 
more freely. 




Figure 15. (Color online) Emergence of a Boundary Breather (BB) following the 
excitation of mode A: = 36 with T = 0.63. (a) time-averaged site kinetic energies (the 
cooling is applied to the first 4 beads); (b) projections of the DB normalized velocity 
field on the NMs as a function of time (non-dissipative dynamics) . 




Figure 16. (Color online) Transient multibreather-likc mode, localized on the protein 
boundary and on the third turn, obtained following the excitation of the fc = 36 NM 
at r = 0.42 (non-dissipative dynamics), (a) Space-time plot of the total site energies, 
(b) Normalized projections of the DB velocity field on the NMs as a function of the 
time. 
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4. Conclusions 

In this paper we liave reported tlie existence of discrete breathers in an off-lattice model 
of proteins with realistic interaction potentials and three-code amino-acid sequence. To 
our knowledge, this is the first time that this class of peculiar vibrational modes are 
found in such a complex and heterogeneous dynamical system. In particular, we have 
shown that, due to the softness of the leading nonlinearity, a necessary condition for 
DB existence is that the linear spectrum exhibits a gap (at least one), as it is indeed 
the case in all proteins. 

For the particular structure analyzed in this work, we have identified three families 
of DBs, each one localized on a different turn of the native fold. The largest fraction of 
the energy of a discrete breather typically resides on a few sites and involves essentially 
harmonic and angular degrees of freedom, while dihedral and Lennard- Jones interactions 
between distant sites are practically not excited. 

We have obtained discrete breathers as continuations of the lower optical band-edge 
normal modes. Quite generically, however, we have shown that members of the three 
different breather families can be excited by feeding energy to any of the normal modes 
lying in the optical band. Remarkably, the excitation of a DB obtained by feeding 
energy to a given normal mode is an extremely efficient process. This means that there 
is an extended range of initial energies that can be fed to the system and immediately 
channeled almost entirely into a breather, which is permanently localized at a very 
specific location. 

Our results prove that DBs have a stabilizing effect on the protein. In fact, a broad 
class of perturbations of the native fold result in the formation of a DB, that is, a 
stable and highly localized vibration, while on the other hand long-range contacts are 
almost devoid of energy. This enhances the structural stability of the protein. To this 
regard, it would be interesting to connect the structure-dynamics relation underlying DB 
formation to the well-conserved patterns of nucleation sites along folding pathways [63] . 

We also discovered peculiar DB-assisted long-range energy transfer phenomena, 
whereby energy is channeled to a distant region of the structure away from the excitation 
site. At variance with ordinary energy exchange between resonant normal modes [M] . 
this is a one-way process, meaning that the transferred energy is never released back to 
the starting location. In particular, this effect may be of importance in the functioning 
of allosteric proteins [HSl El EH] and surely deserves further investigation. 

Our work has highlighted nontrivial correlations between structural and dynamical 
features of protein folds. In a simple /3-barrel, as the one here considered, localization 
occurs on loops. Concerning more complex structures, our results raise the important 
question whether more complex structural selection rules exist driving nonlinear energy 
localization and long-range targeted transfer. 
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Appendix A: Single Site Energy 



Since each potential energy term represents a many-body interaction, we have estimated 
the single-site potential energies V{i) by equally redistributing the energy to each site 
involved in the different terms and then by summing up all this contributions. 

In particular, the harmonic term is a two body potential involving nearby sites, 
therefore the contribution is estimated as follows 



VH{i) 



Y^^2iVH{Tk-i,k)/2 for i = 2,L 
16/(ri,2)/2 for i = l 

Vff(rL_i,L)/2 for i = L 



The angular term is a three body interaction involving three consecutive sites, and 
therefore 

ElltiKi(^.)/3 for i 

\4(^^2)/3 for I 

Va{i) = { EL2VA(^^fe)/3 for I 

T.k=L-2VA{0k)/^ for I 

VA{eL-i)/^ for I 



3,L- 

1 

2 

L- 1 
L 



The dihedral angle <^ is the angle between two nearby planes each containing three 
consecutive sites, the corresponding potential term is therefore a four-body term, and 
the single site contribution can be evaluated as follows 



VDii) 





for 


i = 4,L- 




for 


i = 1 


EL2 Voi^k, Ok, ^fe+i)/4 


for 


i = 2 




for 


i = 3 




for 






for 




Vd{^l-2,9l-2,9l-i)/4. 


for 


i = L 
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The Lennard- Jones term is a two-body interaction involving sites separated by more 
than two bonds along the chain 

VLj{t) = E VLAn,)/2 for \z-j\>2 

3 

The total site potential energies read 

V{i) = Vnii) + VA{i) + VD{i) + VUi). 
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